OPTIMIZED SCHWARZ WAVEFORM RELAXATION AND 
DISCONTINUOUS GALERKIN TIME STEPPING FOR 
HETEROGENEOUS PROBLEMS. 
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Abstract. We design and analyze Schwarz waveform relaxation algorithms for domain decom- 
position of advection-diffusion-reaction problems with strong heterogeneities. These algorithms rely 
on optimized Robin or Ventcell transmission conditions, and can be used with curved interfaces. 
We analyze the semi-discretization in time with discontinuous Galerkin as well. We also show two- 
^s^j , dimensional numerical results using generalized mortar finite elements in space. 
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1. Introduction. In many fields of applications such as reactive transport, far 
field simulations of underground nuclear waste disposal or ocean-atmosphere coupling, 
models have to be coupled in different spatial zones, with very different space and time 
scales and possible complex geometries. For such problems with long time computa- 
' tions, a splitting of the time interval into windows is essential, with the possibility to 

use robust and fast solvers in each time window. 

The Optimized Schwarz Waveform Relaxation (OSWR) method was introduced 
for linear parabolic and hyperbolic problems with constant coefficients in [4]. It was 
analyzed for advection diffusion equations, and applied to non constant advection, 
in |17] , The algorithm computes independently in each subdomain over the whole 
time interval, exchanging space-time boundary data through optimized transmission 
operators. The operators are of Robin or Ventcell type, with coefficients optimizing 
^sO ■ a convergence factor, extending the strategy developed by F. Nataf and coauthors 

O [12] . The optimization problem was analyzed in [5] , [1] . 
• This method potentially applies to different space-time discretization in subdo- 

mains, possibly nonconforming and needs a very small number of iterations to con- 
verge. Numerical evidences of the performance of the method with variable smooth 
coefficients were given in [17\ . An extension to discontinuous coefficients was intro- 
duced in [6 , with asymptotically optimized Robin transmission conditions in some 
k> \ particular cases. 

The discontinuous Galerkin finite element method in time offers many advantages. 
Rigorous analysis can be made for any degree of accuracy and local time-stepping, and 
time steps can be adaptively controlled by a posteriori error analysis, see |20 [ 1X4 } IT6]. 
In a series of presentations in the regular domain decomposition meeting we presented 
the DG-OSWR method, using discontinuous Galerkin for the time discretization of the 
OSWR. In [2], [Sj, we introduced the algorithm in one dimension with discontinuous 
coefficients. In [10], we extended the method to the two dimensional case. For the 
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space discretization, we extended numerically the nonconforming approach in [5] to 
advection-diffusion problems and optimized order 2 transmission conditions, to allow 
for non-matching grids in time and space on the boundary. The space-time projections 
between subdomains were computed with an optimal projection algorithm without 
any additional grid, as in [8j. Two dimensional simulations were presented. In [TT] 
we extended the proof of convergence of the OSWR algorithm to nonoverlapping 
subdomains with curved interfaces. Only sketches of proofs were presented. 

The present paper intends to give a full and self-contained account of the method 
for the advection diffusion reaction equation with non constant coefficients. 

In Section 2, we present the Robin and Ventcell algorithms at the continuous level 
in any dimension, and give in details the new proofs of convergence of the algorithms 
for nonoverlapping subdomains with curved interfaces. 

Then in Section 3, we discretize in time with discontinuous Galerkin, and prove the 
convergence of the semi-discrete algorithms for flat interfaces. Error estimates are 
derived from the classical ones [20J. 

The fully discrete problem is introduced in Section 4, using finite elements. The in- 
terfaces are treated by a new cement approach, extending the method in [8]. Given 
the length of the paper, the numerical analysis will be treated in a forthcoming paper. 
We finally present in Section 5 simulations for two subdomains, with piecewise smooth 
coefficients and a curved interface, for which no error estimates are available. We also 
include an application to the porous media equation. 

Consider the advection-diffusion-reaction equation in il = 

d t u + V • (bu - vVu) + cu = f infix(0,T), (1.1) 
with initial condition 

u(0,x) = u (x) x G 0. (1.2) 

The advection and diffusion coefficients b and v , as well as the reaction coefficient c, 
are piecewise smooth, the problem is parabolic, i.e. v > Vq > a.e. in M. N . 

Theorem 1.1 (Well-posedness and regularity, [13]). Let ft = R N . Suppose 
b G (W^°°(n)) N , v G W 1 ' 00 ^) and c G L°°(fi). If the initial value u is in H 1 ^), 
and the right-hand side f is in L 2 (0, T; L 2 (fl)), then there exists a unique solution u 
of ([TT]) . in jff 1 (0,T;I 3 (n)) n L°°(0, T; H^Cl)) n i 2 (0, T; H 2 (Cl)). 

We consider now a decomposition of fl into nonoverlapping subdomains G 
[1, IJ, as depicted in Figure ITTTl In all cases the boundaries between the subdomains 
are supposed to be hyperplanes at infinity. 




Fig. 1.1. Left: decomposition with possible corners (Robin transmission conditions), right: 
decomposition in bands (Ventcell transmission conditions) 

Problem f| 1 . 1 [) is equivalent to solving / problems in subdomains flj, with trans- 
mission conditions on the interface T^j between two neighboring subdomains fli and 
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flj, given by the jumps [u] = 0, [(z/Vu — b) ■ n s ] = 0. Here rii is the unit exterior 
normal to f2j. As coefficients v and b are possibly discontinuous on the interface, we 
note, for s € j/;(s) = lim e ^.o + K s ~~ The same notation holds for b and c. 

To any z g [1, m], we associate the set A/i of indices of the neighbors of fi^. 

Following [HHi], we propose as preconditioner for (jl.ll 11.21) . the sequence of 
problems 

d t u k l + V ■ (feu? - i*V«*) + au\ = f in iU x (0, T), (1.3a) 

{vid ni - bi ■ n l ) u\ + S itj Ui = (vjd nt - bj ■ n*) u k f x + Siju^ 1 on Fij, j € A/i- 

(1.3b) 

The boundary operators Sij , acting on the part Tij of the boundary of tti shared by 
the boundary of Qj are given by 

Sijip = pijip + qij(d t (p + Vvtj ■ - s.jVr.j^)). (1-4) 

Vr and Vr- are respectively the gradient and divergence operators on T. Pi.j, qi.j, Sij 
are functions in L°°(Ti t j) and r^j is in (L^iTij)) 1 *- 1 . The initial value is that 
of wo in each subdomain. An initial guess (gij) is given on L 2 ((0,T) x Tij) for 
i € G A/i. By convention for the first iterate, the right-hand side in (| 1 . 3[) is 

given by g^j. Under regularity assumptions, solving (11. ip is equivalent to solving 

d t Ui + V • (biUi - v % Vui) + CiUi = / in O, x (0, T), 

(yid ni - bi ■ Bi) Ui + SijUi = (vjdn z - bj ■ B;) Uj + <S,.,". on Tij x (0, T), j e A/i, 

(1.5) 

for i G [1,/] with itj the restriction of u to f2j. 

2. Studying the algorithm for the P.D.E. The first step of the study is to 
give a frame for the definition of the iterates. 

2.1. The local problem. The optimized Schwarz waveform relaxation algo- 
rithm relies on the resolution of the following initial boundary value problem in a 
domain O with boundary T: 

d t w + V- (bw - v^w) + cw = f in O x (0, T), 

vdnW-b-nw + Sw = g on T x (0,T), (2.1) 

w{-, 0) = uq in 0, 

where n is the exterior unit normal to O. The boundary operator S is defined on 
r = dO by 

Sw — p w + q (d t w + Vr • (rw — sVrw)). (2.2) 

The domain O has either form depicted in Figure ITTT1 left for q = 0, right otherwise. 

The functions p, q and s are in L°°(T), and r is in (L°° (T)) N _1 . Either q — 0, and 
the boundary condition is of Robin type, or we suppose q > qo > and the operator 
will be referred to as Order 2 or Ventcell operator. In the latter case, we need the 
spaces 

H' a {0) = {v€H'(0),v\ r eH'(T)}, (2.3) 
which arc defined for s > 1/2, and equipped with the scalar product 

(w, v) H s {0) = (w, v) H s (0) + (qw, v) H s (r) . (2.4) 
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We define the bilinear forms m on and a on H\{0) by 

m(w,v) = (w,v) L 2( ) + (qw,v) L 2( r) , (2.5) 

and 

a(w,v) := I ( — ((b-\7w)v — (b-\7v)w))dx+ j v\7w-Vvdx + j (c + — V • b)wv dx 
Jo 2 Jo Jo 2 

((p — — + ^Vp ■ r)ro» + ^(Vr • (rw)v — Vp ■ (ru)iu) + qsV T w ■ Vrf) dcr, 

(2-6) 

By the Green's formula, we can write a variational formulation of (|2.ip : 

^m(w,u) + a(w,v) = (f,v) L 2 {0 ) + (g,v) L 2 {r) . (2.7) 

The well-posedness is a generalization of results in [Sj [TJ It relies on energy 

estimates and Gronwall's lemma. 

Theorem 2.1. Suppose G W^°°(0), b G (M /1 ' 00 (0)) Ar , c G L 00 (0) ; p G 
L°°(T), q G L°°(T), r € (W 1 ' 00 ^))"- 1 , s € VK 1 ' 00 ^) witfi s > a.e. 

Ifq = 0, iff is inH x {Q,T;L 2 {0)), u is in H 2 {0) and g is in iT x (0, T; £ 2 (r)) n 
L°°(0, T; iJ 1 / 2 (r)) ; satisfying the compatibility condition v d n uo — b ■ nuo + puo = 
g , the subdomain problem (|2.ip has a unique solution w in L°°(0, T; H 2 (0)) D 
W l '°°(0,T;L 2 (O)). 

If Q > qa > a.e., if f is in H 1 ^, T; L 2 (0)), u is in H 2 (0), and g is 
in £r 1 ((0 1 T);.L 2 (:r)) ) problem ((2TTJ /las a umoite solution w in L°°(Q,T;H$(0)) n 
W n > oo (O,T;L 2 (0)) uratfi 5 t u; G L°°(0, T; L 2 (r)). 

Proof. The existence result relies on a Galerkin method like in |18[ fT9] . In the 
sequel, a, j3, ■ ■ ■ denote positive real numbers depending only on the coefficients and 
the geometry. The basic estimate is obtained by multiplying the equation by w and 
integrating by parts in the domain. We set \\w\\ — ||w|Il 2 (C) an( i ll'HIr = ll u 'l!L 2 (r)- 

\ltt m((W ' W ^ + a ( w ' w "> = (/' w )i 2 (o) + (9, w )^(r)- 
With the assumptions on the coefficients, we have 

Case o = 0. 

a(w,w) = j v\Vw\ 2 dx+ j (c + • b)w 2 dx + j ((p - ^-^)w 2 da 
Jo Jo 2 J T 2 

> a(||VH| 2 -0(|M| 2 + Hlr) 

> f I|VHI 2 - 7 IHI 2 , 

the last inequality coming from the trace theorem 

\\w\\ 2 T <C\\Vw\\\\w\\. (2.8) 
We obtain with the Cauchy-Schwarz inequality 

lTt\\M\ 2 + fllVHI 2 < v\\w\\ 2 + S(\\f\\ 2 +\\gf r ), 
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We now have with Gronwall's lemma 

iMi)ii 2 + ? Aiiv^( S )ir-/.s ■ 



2 Jo (2.9) 
2 " T (H|| 2 + 2<5(||/|| + \\g\\ L 2 (0,T;L 2 (r))))- 



We apply (EH to w t : 



IM*)ll a / l|V Wt ( S )|| 2 d S < 

L 2 (UT;L 2 (e>)) + ll<?t|| L 2 (0,T;L 2 (r))))- 

Thanks to the compatibility condition, ui t0 can be estimated, using the equation, by 
IKqII < C(I|woIIh 2 (c>) + 11/(0, Oil) ) and we obtain 

rt 



lM*)lr + 2 y o iiv^(«)irds< 

ff 2 (C>) + (ll/llffi(0,T;L 2 (O)) + NI 

Case <7 > Qo > a.e. 

a(w 7 w) = [ v\\Jw\ 2 dx + [ (c + • b)w 2 dx 
Jo Jo 2 

+ J ((p - + |v r ■ r)«7 2 + g s|V r w| 2 ) do", 
> a(||Vu;|| 2 + ||V r w||r) - 0m(w,w), 

and by the Cauchy-Schwarz inequality 

(f,w)dx+ f (g,w) r da< 1 m(w,w)+S(\\f\\ 2 + \\g\\ 2 ). 
10 Jr 

Collecting these inequalities, we obtain 

±£ t m(w,w) + a(\\Vw\\ 2 + \\V r w\\ 2 ) < (/3 + <y)m{w, w) + S(\\f\\ 2 + \\g\\ 2 ). 
We now integrate in time and use Gronwall's lemma to obtain for any t in (0, T) 



m(w(t), w(t)) + 2a ^ (|| Vw(s)\\ 2 + \\Vrw(s)\\ 2 r )ds < 

e (/3+7)T (m(w , u ) + 25(\\f\\l 2{ ^ T . L 2 {0)) + ||ff||i2 (0)r . i 2 (r)) )). 
We apply (|2~TU|) to w t : 

m{w t (t),w t (t)) + 2a f (|| Vw t (s)|| 2 + \\V r w t {s)\\l)ds < 
Jo 

e ^+-y) T (m(w t0 , w t0 ) + 2S(\\f t \\ 2 L2{0 T . L2{oy) + ||.9i|l|2 (0 ,T;L 2 (r))))- 

We now use the equations at time to estimate m(u; t Q, iy t0 ). From the equation in 
the domain, we deduce that 

IKoll <C(KI| ff 2 ( o) + ll/(0,-)ll), 
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and from the boundary condition that 

IKollr < r]{\\u \\ H 2 [o) + \\g(0, -)||r), 

which gives altogether 

m(w t (t),w t (t)) + 2a f (||Vw 4 (s)|| 2 + \\V r w t (s)\\l)ds < 

Jo (2-11) 

6 e ^ T (\\ u o\\h$(o) + II/IIhMo, T:L 2 (o)) + IIsIIh^o.T;^^))))- 

We can now apply the Galerkin method. When q — 0, we work in i? 1 (0, T; 
H W l '°°(0,T;L 2 (O)) , while if q > q Q > a.e we consider iJ^O, T; H\(0)) H 
W 1>OO (O,T;Z, 2 (0)) n W 1>oo (0,T;£ 2 (r)). This gives a unique solution w. The regu- 
larity H 2 is obtained for g = by the usual regularity results for the Laplace equation 
with Neumann boundary condition, since 

-Aw = -(/ - wt - V • (bw) + Vv ■ Vw) G £°°(0,T; i 2 (C)), 

i" 1 
6>„w = -(b-n-p)w + -q G L~(0,T;£r 1 / 2 (T)). 

In the other case, we have that 
-Aw = -(/ - w t - V ■ (frit;) + • Vw) G £°°(0, T; i 2 (C)), 

vd n - qsArw = i(6 • n - p - q (d t + V r ■ r - (V r • *)V r )))u + i ff G L°°(0, T; i 2 (r)), 
and we conclude like in fTH [19] . □ 



2.2. Convergence analysis for Robin transmission conditions. We sup- 
pose here the coefficients qi to be zero everywhere. Given initial guess (gij) on 
i 2 ((0,T) x Ti.j) for i £ [1,/], j G Mi, the algorithm reduces in each subdomain to 

d t ul + V • (fciuj - i/iVu?) + CittJ = / in a* x (0, T), (2.12a) 

(^A, - 6i • n») +pij = (^-5 n! - 6j ■ Jij) m^ 1 +pij w^" 1 on r ii3 -, j G A/i- 

(2.12b) 

The well-posedness for the boundary value problem in the previous section permits 
to define the sequence of iterates. We now consider the convergence of this sequence. 

Theorem 2.2. For coefficients pij such that Pij + Pj.i > a.e., the sequence 
(uf)fceN of solutions of (|2.12| converges to the solution u of problem (|1.1[) . 

Proof. By linearity, it is sufficient to prove that the sequence of iterates converges 
to zero if / = uq = 0. 

We multiply (|2.12a[) by uf, integrate on f2j, and use the Green's formula. We 
obtain 

5-dl«*(V)llia ( nO + teVu£,Vu£)£» ( n 4 ) + ((<* + • 6,)^,^) L 2 (0i) 



E / (^ 1 uf-^tt, fc )« i k *7 = 0. (2.13) 
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We use now 

(fid^u 1 - - h ■ riiu!- + Pijuf) - (i/id ni Ui - bi ■ niu\ - p %l uf) 2 = 

2(p i , j +p j%i )(u i d ni u k i - *LJ%*)u\ + (p id +p ] , i )(p l , ] - P] , l ~b l -n i )(u l lf. (2.14) 

We replace the boundary term in (|2.13|) . and integrate in time. Since the initial value 
vanishes, we have for any time t, 

\\uUt)\\h ini) +2j\(^utVu*) LHnz) + {(a + iv ■b l )utu k l ) LHnt) ) dr 

+ V / / i (vid ni ui - b t ■ - p %l uf) 2 da dr 

jtC; ■>■■■■ -A /'-,/ • Pi,i 

= 7 (v>id ni Ui -bi-riiUi +p itj Ui) 2 dadT 

jeAr .Jo Jr t}] Pi,j +Pj,i 



E / / (Pi,i - Phi ~ b * ' ni)(u^) 2 da dr. 



Since the coefficients are all bounded, the last term in the right-hand side can be 
handled by the trace theorem (|2.8[) to be canceled with the terms in the left-hand side 
like in the proof of Theorem 12.21 We further insert the transmission condition in the 
right-hand side: 



t 

l«i(*)lli»(nO+"»> ' 



> / W^\\l Hn) dT+J2 f f j (vA^-bi-mu^-pj^u^dadT 

JO jeN -. JO JTij Pi,j +Pj,i 

< V / / \ (v j d ni u^- 1 -b r n i u*- 1 +p i ju'!:- 1 ) 2 d(TdT+C 1 [ 

■ rKr JO Jr< , Pi,j +Pj,i Jo 



\ u i llz, 2 (Q»)- 



We sum on the subdomains, and on the iterations, the boundary terms cancel out 
except the first and last ones, and we obtain for any t G (0, T), 



E E {hHmhi^+vo /Wiii^^) <«(*)+ci E E f 

fceli,x]ie[i,/] v J ° J keli,KJieli,il Jo 

with 



\\ u i lli 2 (a<)) 



(2.15) 



a{t) — ( j (vjd ni u ( * — bj ■ riiU® + pi jii^da dr. 



We now apply Gronwall's lemma and obtain that for any K > 0, 

E E ll«f(*)]li»(n.) < «(T)e^ r , 
fee[i,K] ie[i,i] 

which proves that the sequence u k converges to zero in L 2 ((0, T) x f^) for each i, and 
concludes the proof of the theorem. □ 

Remark 2.3. In the case V ■ b — 0, if Pj,% — Pij —bi -rii = and a > ctQ > 0, 
then C\ = in (|2.15p and we conclude without using Gronwall's lemma. 
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2.3. Order 2 transmission conditions. 

Theorem 2.4. Assume Pij € W 1 ' 00 ^), p uj + Pj,i > a.e., q it j = q > 0, 
h e {W 1 ' 00 ^)) 1 *, Vl e W^iSi), r itj e (W rl - 00 (^)) A ' Jl , wtti»r w = r jti on r h3 , 

Sij € T / F 1 '°°(r2i), Si.j > with Si.j = Sjj on Tij, and the domain is cut in bands as 
in Figure right. Then, the algorithm ()1.3j) converges in each subdomain to the 
solution u of problem (jl.lj) . 

Proof. We first need some results in differential geometry. For any i € [1,/], For 
every j € Mi, the normal vector rii can be extended in a neighbourhood of T^j in O,; 
as a smooth function fit with length one. Let ipi t j € C°°(f2j), such that = 1 in 
a neighbourhood of r,-^-, = in a neighbourhood of T^k for k £ Mi, k ^ j and 
V'ij' > on fii. We can assume that Hi is defined on a neighbourhood of the 
support of ipi,j- We extend the tangential gradient and divergence operators in the 
support of ipij by: 

Vr^V : = v ¥> ~ V r ,. 3 ■¥> := V • (<p- {<p- n^fn). 

It is easy to see that (Vr ; ^v)|r< ■ = Vr< j<p, (Vr 4 , • v)|r; ■ = Vr/j ^ ■ and for <£> and 
X with support in supp(ipij), we have 

/ (V rili -<p)xdx = - [ tp-Vr^xdx. (2.16) 

Now we prove Theorem 12.41 We consider the algorithm () 1 . 3[) on the error, so we 
suppose / = w = 0. We set \\(p\\ z = |M|z,2(o, .), ||M|| f = ||^77~V<^||f, IMkoo = 

IMU°°(n<)j IMkl>°° = H^ll 1^1.00(0,) and = yJ ?M±PM , 

The proof is based on energy estimates containing the term 

/ / {ndnt u k - h ■ n t u k + Siju 1 *) 2 da dr, 

JO Jrij 

and that we derive by multiplying successively the first equation of (|1.3p by the terms 
pful, d t u k , Vr 4i< • and -Vr 4l< ■ V r « d uj). 

We multiply the first equation of f| 1 . 3[) by 0fu k , integrate on (0,t) x and 
integrate by parts in space, 

\mu k S)fi+ flWiuHr,-) III? dr- f [ ■ V fr)^) 2 dx dr 

L Jo Jo Jn z 

+ 11 (c 8 + ^V ■b i )l3j{u k l ) 2 dxdT~ f f Ui\V0i\ 2 (u^ 2 dxdT 
Jo Jo.i * Jo Jili 

flA"* 9 *** - ^u^u^dadr = 0. (2.17) 

We multiply the first equation of (|1.3p by dtu k , integrate on (0,t) x f2j and then 
integrate by parts in space, 

I \\d t u\\\ 2 dT+\\\\u k S)\\\ 2 + I [ (c l u k +V-(b l u k ))d t u k dxd T - f f vA t u k d t u k dadT = Q. 
Jo * Jo Jili Jo JVi j 

(2.18) 
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We multiply the first equation of (I1.3P by Vr ; s ■ (ipf j r i,j u i) integrate on (0, t) x 
and integrate by parts in space: 

dtu* Vr itj ■ (ipljTiju^) dxdr+ [ [ V • (b t u k ) Vr^ • {^l j r l , 3 u k ) dx dr 

'0 JUi Jo JUi 

t 

CittiVr 4ij • ($1^11$) dx dr 

'0 Jfi 



-[ I v i VUi-VVr i , j ii>ljri,ju!?)dxdT- f [ Vid ni u k Vr itj iri,jU k ) dadr = 0. 

JO JUi JO JTi j 

(2.19) 



We observe that 



v l Vu k ■ VVr^. • (ipl 3 r lt3 u k )dxdT 



= [ [ ^u k -V(Vr t y(^l :j r l:3 )u k )dxdT+ f f ^Vi^-V^Vij-Vr^ u k ) dxdr, 
Jo Jut Jo Jrii 

(2.20) 

with 

/ / i/iVu? • V(^r 4J • Vr 4J u?) da; dr 
Jo Jn ; 

>-7 ril^VTOVVr^ttfllfdr-C / / (||V«?||? + ||«?||?) dxdr. (2.21) 
4 Jo Jo JUi 

Replacing ([2~2~T|l in (|2~20ll and then (f2~20|) in ((2~T9|) . we obtain 

^9t«i V r<ij ■ {ipljrijUi) + V ■ (fcjtii ) Vr Ii3 • (ipi,jri,jU k ) + au k V r<J • (^-r^u*)) dxdr 



JO JT; , 

t 

fe||2 , no „Jc||2\ 



<C / (||V^V^|| 2 + ||A U f||^)dxdr. 

JO Jfi; 

(2.22) 

Now we multiply the first equation of f) 1 .3[) by — Vr; 3 - • (ipf jSij Vr; 3 . w^) integrate on 
(0, t) x £!j and integrate by parts in space: 



Jfi; 

+ / / V'i.jSi.iVr^CciU*) • Vr id u k dxdr- / / fiVu* • V(Vr 4iJ • {i'ijSij Vr 4 ^ dsdr 

JO Jf!; JO Jf2; 

+ / / ^d„A fc Vr i3 - •(si,jV r<3 ^)d CT dr = 0. (2.23) 
Jo Jr, , 
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We have, 

-If vtfu\ ■ V(V r<J ■ WijSij V riJ . u?)) dx dr 

>IJ Q UijV^U V Vr 4>j «?)||? dr - d jf ||V«?||? dr. (2.24) 
Replacing ([2~2"4"]l in (|2~2"5)) leads to 

^ll^V^JVr^ «?(*)||? + \J UijV^IJ VVr l>3 u?)||?dr 

+ 11 ^l j s i ,jC l \Wr lJ u'l\ 2 dxdT+ I / i/^tt* V r , j ■ (s itj V r , 3 u^) da dr 
Jo JQi ' ' Jo Jr itj 



< 



I I V • (fcu* ) Vr itj ■ ($j8i d V r , , u\) dxdr + C f ||v^V«? ||? dr. (2.25) 

JO J S~2; JO 



Multiplying (|2~T5)l . (|2~2"2")l and (j2T2"5)) by q, and adding the three equations with (|2TT7|) . 
we get 

5 (||A«?(*)II? + 9 III u?(t) III? +g||^iJV^7 V riJ «?(t)||?) + j |||ft«?(r, •) III? dr 

Pf(v t d n ^ - bj ^u'l) uj dadr 
'o Jr 4 ,i ^ 

/• 

i/i^uj (ftuj + Vr,,, • (r^uj) - Vp M • (s^ V r ^ u, fe )) t&rdT 

o Jr*^ 

< «(|||6i|N,i,oo + IKilki,*,) I ||ti*||< llftt^Hi rfr 

+ «(||6i||i,oo + ||rij|| 1 ,oo) / ||V«?||i||^||(dr 
Jo 

+ |l|6ill*,oo y ||V«J||i||Vr 4iJ •(^VijVr 4J ttf)|| i dT 
1 f* 

- sCglNki.oo + ||c<IU,oo) / IK II* ||v r<iJ • {i>ijSi,j v r ,, dr 





+ I (|Nki,°c + ||ci|koo)||«?(t)||? + C ( f ||A«?||?dr + g / " HV^V^IIi ^ 



JO 



We bound the right-hand side by 



+ |(l|6|||i,i,oo + ||cilkoo)K(t)||f + ^( / ||ftuf|kdr + g/ H^VufH^r 
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We simplify the terms which appear on both sides, and obtain 

\ (||A«?(*)||? + 1 III u?(«) III? +?II^V5«Vr 1J «?(t)||?) + jT* |||ft«*(r, •) |||? dr 
2 Jo 



'o ^r 4 ,i 2 
fi9 ni u^ (<9 t ttf + Vr,,, • (rij-uj) - Vr l>3 • (s fjj V r<>j u%)) dadr 

<cQ \\fcu*\\UT + qJ \\V^vZ\\Ut) ■ (2.26) 
Recalling that s^j — Sjj on Tij and = rji on T^j, we use now: 



= 4 (^P'lj{vid ni u k i - bl J** Ui)ui + qfid ni Ui(d t Ui + Vr M ■ (r»,jU< ) - Vr M ■ (s;,j Vr M «>))j 

+ 2<?(p iJ - p Ji( - 26, • niXftuJ + V r «j ■ (r w u?) - V r , :] ■ (s lJ Vr I>3 

+ (pij + Pj,i)(pi,j - Pj,i - 6; • n i )(M l fc ) 2 . (2.27) 

Replacing (j2~2T)) into we obtain 

5 + 9 III «?(*) III? +9ll^-^JVr 4J «?(<)!!?) + jT \muUr, •) HI? 

+ |/ \\ d t u i\\i dT + \ ( ( {v l d n ,u k ; -bi-tiiU* - Sj,iuf) 2 dadr 

8 Jo 4 Jo Jr^ 

+ / / iPi,j+Pj,i)(-Pij+Pj,i + b i-ni)(uf) 2 dadT+ ^(||6i||i,l,oo + l|c<ll»,oo)]|«Ai (*)||f 

Jo Jrij ^ 

+ tfQf ||A^||^r + Q^ Hv^Vufll^dr) . (2.28) 
In order to estimate the fourth term in the right-hand side of (|2.28j) , we observe that 

n(-Pi,j + Pj,i + 26, ■ n^u* d t u\ dadr = ^ / (-p M + + 26, ■ n^u^t) 2 da. 
-'i,3 JTiJ 

By the trace theorem in the right-hand side, we write: 



10 JTi 



(~Pi,j +Pj,i + 2bi ■ Ui)u^ d t u- dadr < (7||<(t)|| i || v ^'V<(t)||i, 



12 



L. Halpern, C. Japhet, J. Szeftel 



and , 

ft 

,k\„,k ^ n I I 1 1 n . k 1 1 2 i / / 1 1 „ , k 1 1 2 



<(*)ll? = 2/ / (9 t <)<<2(/ WdtufWi) (/ ||« > (2-29) 

JO JQi 



we obtain 



| / / (-Pt,3 + + 2&i • 9 t M 4 fe da dr 

2 Jo Jr, ,• 



< | jf ||ft«*||?dr + I III «?(t) |||? +C ^ ||A^||?dr) ■ (2-30) 
Moreover, integrating by parts and using the trace theorem, we have 
_ o/ / v r i3 • (si,j ;Vr 4 , Ui)(-Pi,j +Pj,i + %bi ■ n^u* da dr 

< ^ y II^V^J VVr i!3 . uf ||? + C^' ||V r<J u?||? dr + ^ ||A«?||? dr). (2.31) 
Using (|2.29j) . we estimate the third term in the right-hand side of (|2.28[) by 

|(HMlM,=o + IMIi,oo)||u?(t)||? < § j ll^i II? dr + C j ||ft^||?dr. (2.32) 

Replacing f|2.31[) . (|2.30[) and (|2.32[) in (|2.28[) . then using the transmission conditions, 
we have: 

\ (l|Au?(*)ll? + f III «?(*) Ill, +?H^,iVSiJVr 4iJ u*(t)||? 

' |||/3^ fe (r, •) III? dr + | y' ||9 t ^ fc ||? dr + ± J* HVwVTO V V r ,. 3 ^|| 2 dr 



1 '* 



4 7o 



9 ni u^ — bi -rii u\ — Sj^uf) da dr 



1 



- - / / (vjd ni vh l -b -n l u k 1 1 +S itj Uj : ) dcrdr 



+ c(y t ||^ fe || 2 dr+| " 



y iiv^v^ fc n?dT 



We now sum up over the interfaces j G A/i, then over the subdomains for 1 < i < J, 
and on the iterations for 1 < k < K , the boundary terms cancel out, and we obtain 
for any t G (0,T), 



E E (ll/8iJ«i(*)lli+?llv^Vu?(t)||? + i^y 



||V(A,,-^)||? dr 



fee[i,if]ie[i,/] 



<a(i)+C E E (y HA,i^||?dr + g y ||V^V^||?dr), (2.33) 



fce[i,K]*e[i,il 

with 



m = ~ e E / / fas*"? - ^ • ^ m ° + s i,3 u °o) 2 d ° dr - ( 2 - 34 ) 



We conclude with Gronwall's lemma as before. □ 
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3. The discontinous Galerkin time stepping for the Schwarz waveform 
relaxation algorithm. In the following sections, in order to simplify the analysis, 
we suppose that c + • b > ao > a.e. in tt. 

3.1. Time discretization of the local problem: discontinuous Galerkin 
method. We suppose that the coefficients are restricted to p - ^ + | Vr ■ r > 
a.e. on T, q > a.e. and s > a.e.. This implies that the bilinear form a defined in 
(|2.6p is positive definite on when q = 0, and positive definite on H\(Q) when 

Q > qo > a.e. 

We recall the time-discontinuous Galerkin method, as presented in [14]. We are 
given a decomposition T of the time interval (0, T), /„ = (t n ,t n _|_x], for < n < N, 
the mesh size is k n = t n+ i — t n . For B a Banach space and / an interval of R, define 
for any integer d > 

d 

P d (S) = {y> : 7 -> 6, = W**' ^ e 

T) - {</? : I -> 6, ^ G P rf (B), 1 < n < N}. 

Let B = H{(0) if q > 0, 6 = .ff^C) if g = 0. We define an approximation U 
of u, polynomial of degree lower than d on every subinterval /„. For every point 
t n , we define U{t~) = limt_).t n _o ?7(i), and note U(t^) = lim t ^,.t n +oU(t). The time 
discretization of (I2.7P leads to searching ?7 € Vd{B,T) such that 

' Cf(0, •)=«<>, 
We V d (B,T): J ( m(J7, V) + a(t/, V)) A 

+m(f/(t+ •) - U(t~, ■), V(t+, •)) = J i L(V) dt, 

with L(V) = (/, V^) L 2( ) + (g, V) L 2(i>\. Since I„ is closed at i„+i, J7(t^ +1 ) is the value 
of U at t n+ i. Due to the discontinuous nature of the test and trial spaces, the method 
is an implicit time stepping scheme, and U G Pd{B, T) is obtained recursively on each 
subinterval, which makes it very flexible. 

Theorem 3.1. If p — ^ + §Vr • r > a.e. on T, q > a.e. and s > a.e., 
equation (|3.ip defines a unique solution. 

Proof. The result relies on the fact that the bilinear form a is definite positive. 
It is is most easily seen by using a basis of Legendre polynomials. U G Vd(H\(®),T) 
is obtained recursively on each subinterval. We introduce the Legendre polynomials 
L n , orthogonal basis in L 2 (— 1,1), with L n (l) = 1. L n has the parity of n, hence 
L n (— 1) = (—1)"- A basis of orthogonal polynomial on /„ is given by L n ,k(t) = 
L k (i^(t - t " +1 2 +t " )). Choose V(t,x) = L n j(t)$j(x) in (JOJl with G #f(ft), and 

expand [/ on 7„ as U(t,x) — Ylk=o^k{x)L n ^{t). Suppose [/ to be given on (0,t n ]. 
In order to determine £/ on /„, we must solve the system: for any <frj G H\{£l), 

/ f /'.„ /./.„ ,//'•:/ /,. <!>, : + L n!k L n ja(Uk,$j) J dt 

d 

+ 5jL n , fc (t+)L nj -(t+)m(U' fc ,$ i )= / LnjL^dt. 

k=0 Jl " 
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It is an implicit scheme. We calculate the coefficients 



/ Ln fcLn j — 8fcj || L n j || j 
Jin 



i£k<j 

1 - (-l) fe +J if k > j, 

'(_l)*+3 iik<j 

1 if k > j, 



L n ,kL n .j + Ln,k(tn)L7i,j(t^) 



which leads to 



\\L n ,i \\ 2 a(Uj,*i) + mQJi^i) + ^(-l) fc+J m(t/ fc , *,■) + £ m(t/ & , *,) = /" L rej L(%) dt. 

It is a square system of partial differential equations, of the type coercive + compact. 
By the Fredholm alternative, we only need to prove uniqueness. Choose now $>j = Uj, 
and obtain 

^\\L n J\ 2 a(U jl U j ) + ^m(U j ,U j ) + 2'£ i m(U k ,U 1 )=0, 

j j j k>j 

k+j even 

and since a is positive definite, we deduce that U = 0. □ 

We will make use of the following remark ([16]). We introduce the Gauss- Radau 
points, (0 < n, . . . , Td+i = 1), defined such that the quadrature formula 

d+l 



Jo „•_•, 



is exact in P2d, and the interpolation operator I n on [t n ,t n _|_i] at points (t n ,t n + 
rife„, . . . ,t n + T d +\k n ). For any \ G Pd, X = InX G Pd+L 

Let I : Vd(B,T) — > Vd+i(B,T) be the operator whose restriction to each subin- 
terval is X n and satisfies IU{b+) = U(t~). By using the Gauss-Radau formula, which 
is exact in P2d, we have for all ipij G Pd 

J I> ~dt~ ^ * ~ Jj Tt ^ dt = ~ 
As a consequence, we have a very useful inequality: 

/(Z 1 
— (ltp id )ip i:j dt > - 

Also, equation (|3.1j) can be rewritten as 



(llPi^TPijdt > -[^Atn+l) 2 ~ TpiAtnf}- (3-2) 



(m(^,V) + a(U,V))dt= [ L(V)dt, (3.3) 
dt J In 

or in the strong formulation: 

d t {2U) + V -(bU- vVU) +cU = Pf, in Q x (0, T), 
(vd n -b-n)U+pU + q{d t (XU) + V r ■ (rE7 - sV r f7)) = P.g on T x (0,T) 



(3.4) 
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Here P is the projection L 2 in each subinterval of T on P^. 

Theorem 3.2 (Thomee, [20] )• Let U be the solution of (|3.ip and u the solution 
of (|2.ip . Under the assumptions of Theorem \2.1[ the estimate holds 

\\u - tf|U- ( j B ,i»( n) ) < Ck d+1 \\d? +1 u\\ LHQ ^ Hi{n)) , (3.5) 

with k — maxo< n <Ar k n . 

3.2. The discrete in time optimized Schwarz waveform relaxation al- 
gorithm with different subdomains grids. In this part we present and analyse 
the discrete non conforming in time optimized Schwarz waveform relaxation method. 

The time partition in subdomain Qj, is 71, with A^ + l intervals 7^, and mesh size 
k\. In view of formulation (|3.4p . we define interpolation operators T l and projection 
operators P l in each subdomain, i.e. P % is the projection L 2 in each subinterval of % 
onPd, and we solve 

d t (TUf) + V ■ (hU? - ViVUf) + Cl L/* = P l f in n t x (0, T), (3.6a) 

( Vl d n% - bi ■ m) C7? + SijU? = P^ivjdn, - bj ■ m) XJ]~ X + s^u]- 1 ) on r id ,j e N t . 

(3.6b) 

Here the operators are different on either part of the "equal" sign: 

S itj U = PiJ U + %j (d t {TU) + Vr 4J • {njU ~ 8ijVr td U)) ,„ ?) 

S id U = jhj U + q itj {8t(pU) + V r „ • (r itj U - a id Vr tJ U)). ' 

Formally, the sequence of problems (|3.6p converges to the solution of 

d t (TU t ) + V • {hUi - ViWUi) + c t Ui = P l f in tt t x (0, T), (3.8a) 

(vid ni - b % ■ rii) Ui + SijUi = P % ((vjd ni - bj ■ n<) Uj + SijUj) on T iyj ,j e M- 

(3.8b) 

We present the analysis first with Robin transmission conditions (e.g. qij — 0) and 
general decomposition, and then with order 2 transmission conditions and decompo- 
sition in strips. 

3.2.1. The Robin case. We consider here a general decomposition of the do- 
main, possibly with corners. We solve (|3.6p with q^j — 0, i.e. SijU = Si.jU — Pij U. 

Theorem 3.3. Assume qij = 0, Pjj — Pij — bi -n,; = 0, Pij — ^p- > 0. Problem 
(|3.8p has a unique solution (£7j)jej , and Ui is the limit of the iterates of algorithm 
(EH)- 

Proof. We first write energy estimates on (|3.6p for / = and uq = 0. We start 
like in the proof of Theorem 12.21 We multiply (|3.6a|) by U-*, integrate on f2j, then 
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integrate on the interval (t l n ,t l n+1 ) and use (|3.2I) and (|2.14l) : 

ll^i (*n+l)lll, 2 (fi«) — (*^)||i2(Q.) 

+ 2 [ ((viVU t k ,VU?) LHQi) + ((a + l -V ■ bi)Ut Ufivw) dT 



+ V / / (vid^u 1 ; - h ■ riiUf - p^Uffda dr 

< V / / (^d^uf - bi ■ n.Jjf + PijUffda dr 

+ V / / {Pj,i-Pi,j ~b l -n l )(U k ) 2 dadr. 

We can not use Gronwall's Lemma like in the continuous case, due to the presence of 
the global in time projection operator V 3 in the transmission conditions. Therefore 
we have to assume that pj^ — pi.j — bi -Jlj = everywhere, which cancels the last term. 
We sum up over the time intervals, using the fact that the errors vanish at time 0: 

l|tf*CT)lli»(no + 2mim>o,a ) / \\Uf \\ 2 m{sii) dr 



V / [ i {vid ni vll - b t -riiUf - Pj, t U k ) 2 da dr 

< V / f {vid ni u\ -bi-niU? + Pi , j U k ) 2 dcrdT. 

r \* JTij Pi,j + Pj,i 



We now insert the transmission conditions 

rT 



||^(T)||| 2(Oi) +2mhn> ,a ) / || U? \\ 2 HHni) dr 

Jo 

+ V / [ 1 (vid ni u k -bi-niUf - PjA Uf) 2 d<rdT 

^Jj-.JO J\ . /'-,; • I',-. 

* E T / ('"(''^' ; : t>mj : 1 • ))\hd-. 

j eJ y.Jo Jr id Pi,j +Pj,i \ J 
We use the fact that the projection operator is a contraction to obtain: 

\\U k (T)\\l 2m +2rmn(^a () ) f \\U?\\ 2 Him dr 

Jo 

+ V / / 1 (vid ni Ui -bi-mllf -pj^fdadr 

je sf.Jo Jv /'• .. • I'j - 
We sum up over the subdomains, we define the boundary term 

B " = E E f I 7 - h ■ mUf - p 3 Mfd<y dr, 
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we obtain 

E (\\Ut(T)\\h m +2mm(^a ) f Wll^n.) *0 + B k < B^ 1 . (3.9) 

We first apply this inequality to prove the first part of the Theorem. (|3.8p is a 
square discrete system, and proving well-posedness is equivalent to proving unique- 
ness. Dropping the superscript in (|3.9[) gives the result. As for the convergence, 
we proceed as in the continuous case by summing p. 91) over the iterates to obtain 

that _E ie p,i] H^Cniliacn,) and Eiep.i] Jo ll^lljr^n,) dl ~ tend to zer0 as k tend to 
infinity. □ 

3.2.2. The Order 2 case. We restrict ourselves to a splitting of the domain 
into strips with parallel planar interfaces. 

Theorem 3.4. We assume that pi^j = p > 0, qij = q > 0, s^j = s > 0, 6j = 
and Tij — 0. Problem (|3.8jl /ias a unique solution (Ui)i£i , and Ui is the limit of the 
iterates of algorithm ()3.6j) . 

Proof. We consider the algorithm (|3.6p on the error, so we suppose / = uq = 0. 
As in the continuous case, the proof is based on energy estimates containing the term 

J J {vidn^+SijUff dadr, 

and that we derive by multiplying successively the first equation of (|3.6[) by the terms 
Ul d t (l%*), and -A r<iJ Uf. We set ||M||f = Wy/viV^^ + ||V^||| W We 
multiply the first equation of (|3.6II by f/f, we integrate on P n x fii then integrate by 
parts in space and use ()3.2[) : 

511^(^1)11?+/ Illttffo Widr-f / z/A ( ^f^T<i||c/f(i-)|^ 

(3.10) 

We multiply the first equation of (|3.6p by d t {l l U^), integrate on P n x fij and then 
integrate by parts in space and use (|3.2[) : 

5111^+1)111?+^ WdtiVUDWjdT-J^ f Vl d n ^d t {VU^)dadT< l -\\\U^t-)\\\i. 

(3.11) 

Now we multiply the first equation of (|3.6I) by — Arv , f/f integrate on 7^ x f2, and 
integrate by parts in space and use (|3.2p : 

5l|Vr 1 , J C/ffc +1 )||f + ^ |||Vr 4J I^(i,.)|||? 

+ "i^^Ui^dT^illVr^^t-JH?, (3.12) 

where we have used the fact that Arv is a constant coefficient operator. Let 
E n {Ut) = §||C/f(«)-)||? + § III E/?((0") III? +^||V r4J t^(«)-)||?. 
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Multiplying (|3.10p by p, (|3 . by q and (13. 12)) by sq, and adding the three equations 
with (I33U1) . we get 



E n+1 (U?)+ ^[P III -)|IL 2 +q\\d t (T n U t k )\\* +sq\\\Vr i , j U*(t,-W i }dt 

~ E I [ Vid "< U ^ 5 « U i dx2 dt < EU ( U ; 

It can be rewritten as 

lfTTk\ , / nurrku \ Ill2 , „ll a /<ri rrfcM|2 , „„ in T7 rrfe/'i mm 2 - 



We now sum in time for < n < N, and use the transmission condition. Since 
E°(Uf) = 0, we obtain 

E N+1 (U?)+ f T [\\\U?{t r )\\\t +q\\d t (T n Uf)\\i +*g|||Vr 4J ^(i,-)|||?]d* 
Jo 

+ lE/ T / ("AHvt-SijU?)*dt<\ f I (P^-Vjdnrf-'+SijUf-^fdt. 

We sum up over the subdomains and use the fact that the projection is a contraction. 
Since we are in the case where pij = p, qij = q, = and Sij = s, we have Sij — Sji- 
Thus, we can sum up over the iterates, the boundary terms cancel out, and we obtain 

|[/^,.)||| 2 + q \\d t (V n U*)\\ 2 l +*g|||Vr 4 X(V)lll?]<&) 



f;E(V + W)+/ 

k=i i=i V Jo 

+ z e r / ^< - s « u ?? < i e c i - saw 



We conclude as in the proof of Theorem 13.41 □ 

We now state the error estimate in the Robin case. 

3.3. Error estimates in the Robin case. Theorem 3.5. If V • b = 0, 

Pi,j — hi 2 "' = Pj,i — hj 2 " J = p > 0, and qi j = 0, the error between u and the solution 
Ui of (|3.8[) is estimated by: 

I 

<Ck^ d+ ^dt + \\\l^ T , m{n)y (3.13) 

i=l 

Proof. We introduce the projection operator Pr as 
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We define Wj = P i (wloj, 6j = E/j — Wi and pi = Wi — u^. Classical projection 
estimates (|20j) yield the estimate on pf. 

i 



L 2 (0,T;L 2 (fi))- 



(3.14) 



Thus, since f7, — u\n t = Oj + Pi, it suffices to prove estimate (|3.13[) for 6^. Now, 
thanks to the equations on u and Ui, and the identity 4j^P^ = T >% 4z, Oj satisfies: 

9 t (I l 0i) + 6 • V9, - uA@i + cQ l = b ■ Vpi + vAp t - cp t 
+ (1- P l ){d t u- /) in SU x (0,T), 
[yi d n , -bi-ni)<di+ pi j Qi = P l ( (i^j 9 ni - b ■ n l ) Qj + pij Qj ) 

-(i - p i ){{v 3 d ni - bj ■ m)Wj + p itj Wj) on rv, x (0, T), j G M- 

Multiply the first equation of (|3.14[) by Oj, integrate on (i^j^+x) x Oj, using p. 21) 
and integrate by parts in space. Terminate with Cauchy Schwarz inequality: 

|ll©<((*«+i)~)ll< + jf llie<(*> ■) III? dt ~ lJ r (yi®ni Q i - ®i dx * dt 

Rewriting the boundary integral using (|2.14| , we obtain 

Jl|0i(«+i)")lli+/ ll|e»(V)|||?*+^- J2 f f {vid ni Q i -b v n i Q i -p j ,Sifdx 2 dt 
1 J i' n 4p jeK Jl » Jr *i 

^ r E I I (^9 ni e i -6 i -n i e i +p iJ e i ) 2 dx 2 +i||e i («)-)||?+c / \\Pi(t,-)f HHQi) dt. 

Using the transmission condition in Q3.14J) together with the fact that P l and 1 — P % 
are orthogonal to each other and have norm 1, we get by a trace theorem 



1 ,2 



dl0*((4+i)-)l 



?+/ |||6 I (t,-)||L 2 d<+^ V / / (i/i^ei-fti-ni^-p^eO'dcadi 



^ i E / / (^^% + |l|0z((4)~)ll? 



+c y^ iipi(*.oiiH>{n 4) ^+c y j4 Ha -^(ttk)^. oii^cno ( 3 - 15 ) 

Classical error estimates [5U] imply: 

IIM*. Oll^cno di + l T \\(l-n(u\n,)(t, Oilmen,) < Cfc 2(d+1) ||9 t d+1 W ||i 2(0 , T;ff2(Oi)) . 

Summing ()3.15p in i and n, and using the previous equation yields (|3.13|) . □ 

4. Space-time nonconforming algorithm. In this section we describe the 
implementation of algorithm (13. 6p . especially in the cases d = and d = 1. We start 
from the semi-disrete in time scheme and use finite elements for the space discretiza- 
tion in each subdomain. In order to permit non-matching grids in time and space on 
the boundary, we extend the nonconforming approach in [8] . 

We describe first the implementation of algorithm (|3.6[) at the semi-discrete in 
time level, and then at the space-time discret level. 
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4.1. Time discretization. We recall the subdomain scheme in time, and give it 
in details for d = and d = 1. Then we describe the computation of the transmission 
conditions in algorithm (|3.6p . 



4.1.1. Interior scheme. We consider the subdomain problem in the algorithm 
(GEU) at iteration kmO = 0;. Let B 4 = /^(fi;) if q > 0, B< = H 1 ^) if q = 0. We 
set U = € Vd{Bi,7i), and we omit the subscript i for the local time scheme to 
simplify the notations : 



d t (!U) + V • (bU - uVU) +cU = Pf in O x (0, T) . . 

(i/<9„ -6-n)f7+jjf7 + g(9t(X?7) + V r • (rU - sV r U)) = Pg on T x (0,T)/ > 



Case d = 

In the case d = 0, the approximating functions are piecewise constant in time, then 
U(t) = U n+1 = U% ml n , wehavel n U = U n + l -=^{U n+1 -U n ), P£ = £ ((-,s)ds 
and the method reduces to the modifed backward Euler method 



Tjn+l jjn i f 

V ■ {bU n+1 - v^U n+1 ) + cU n+1 = — + — / /(■. s ) ds in ° 

tin ti n J J n 

Tjn+X 

P U n+1 + q(—— + X7 r -(rU 

tin 

U n If 
q T~ + T / d(-,s)ds on T. 

tin tin J i 



Tjn+X 

{vd n -b- n)U n+1 + v U n+1 + q(-r— + V r ■ (rU n+1 - sV r U n+1 )) (4.2) 

tin 



Case d = 1 

In that case, for piecewise linear functions of f, using a basis of Legendre poly- 
nomials we may write, 17(f) = U^ l+1 + 2 *~^+ 1/2 t/ 1 " +1 , on with t n+1/2 = , 
U n = U(t n ), and we have on P n : 



l n U = - U{ 1+1 - U n ) + (UZ +1 + [/f +1 - U n )( - tw+1/2 ) 

+3(-U£ +1 + U{ 1+1 + ^")( * ~*" +1/2 ) 2 , 

k n 



and P£ = £ + 2 " +1/2 gi with 



(4.3) 
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Thus, we obtain for the determination of Uq and U™ the system 

1 TT n 1 f 

tt(U 7 o +1 + U[ l+1 ) + V ■ (6f/ " +1 - ^VC/ " +1 ) + cf/ " +1 = + — / /(., s ) ds 

|-(-[/ " +1 + C/ r +1 ) + V ' ( bU ? +1 - ^VC/ 1 n+1 ) + cU{ l+1 
3U n 6 f s-t n+1/2 



J I„ ™n 



h h 

(vd n -b- n)U» +1 + P U r l+1 + q(^(K l+1 + Up 1 ) + V r ■ (rUp 1 - sW r Up v )) 

fan 

1 f , * J 

= Q-j— + T~ / 9{-,s)ds 

(ud n - b-n)U? +1 + P U[ l+1 + g(^-(-t/ n+1 + Up 1 ) + V r • (rUp 1 - sVrUp 1 )) 

3(7" 6" f s - t n+1/2 
= -q-r h — / g(-,s)ds on T. 

(4.4) 

Multiplying the first equation of (|4.2I) by w e 8, (resp. the first equation of (|4.4I) 
by w G ^ and the second equation of (|4.4I) by u> € Si), integrating by parts on O, 
and using the boundary conditions, the variational formulation is: 

Case d — (Variational formulation) 

m(U n+1 ,v) + k n a(U n+1 ,v) = 

m{U n ,v) + (f(;s),v)ds + J^(g(;s),v)rds, V«e6 t . (4.5) 

Case d = 1 (Variational formulation) 

m(L/o n+1 ,u) + k n a(Up\v) + m(Up\v) 

= m(U n ,v) + (f(-,s),v)ds + J (g(;s),v) r ds, 
-m(Up 1 ,v) + m(Up l ,w) + k n a(Up 1 ,w) 

= -m(U n ,v)+f 2{S ~ tn+1/2 \ f(; S ),w)ds 

+ / 2{S ~* n+1/2 \ g(;s),w) r ds, Vv^weB,. (4.6) 



Remark 4.1. Equations (|4.5p and (|4.6p can oe derived directly from (|3.ip . How- 
ever we will need formulas (|4.2p and (|4.4p in t/ie space nonconforming case. 

We now discuss the computation of the right-hand side on the interface Tij x (0, T) 
for j e A/j in the algorithm (|3 . 6[) . 



4.1.2. Transmission terms. Let (#* j) be a given initial guess in Vd(L 2 (Tij), %), 
for 1 < i < J, j 6 Mi. Then, at iteration k > 1, we solve the subdomain problem in 

dtipU*) + V • (hUf - i/iVUf ) + c, = P7 in ^ x (0, T), (4.7a) 
(i/<£U, - h ■ tii) Uf + SijU? = gl v j G M on T 4J x (0, T), (4.7b) 
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The function g*. is denned for k > 2 by 



9i,i = I "til- (4-8) 

with k > 2, defined by 
We remark that, for k > 2, 

&i = -gfc 1 + s J , i u*- 1 +s itj u!i- 1 

= -9 -J 1 + (Pi,j + Pj,i)U^ 1 + (qi d + qj.MFUj 'i 

Once <?j? j is computed from Uj 1 , we obtain from (|4.8p as follows : we introduce 
the basis functions (ip l n a )o< a <d of polynomial of degree lower than d on subinterval 
4, then 



(9i,j)\Ii, = i pl 9j t i)\I* n = ^ Q^n.q 

with GJ; fc Q € L 2 (rjj) solution of the system 

J2 G n% [ rfi,a<Pn,0ds= [ #^,0, /3 G {0, d}. 



Thus, the computation of gfj on each I % n needs the computation of terms in the form 

9j,i<Pn,f> ds ' ( 4 - 9 ) 



for /3 G {0,...,d}. Recall that is defined on T tJ x I and G V d (L 2 (T ifj ),Tj). 
Thus, we first write the integral in (|4.9p as an integral over / : let <£>^ a be the function 
defined on /, equal to ip l n a on I l n and equal to zero on T\l % n . Then 

%k ,J j _ / ~fc ,f.i 



= / fl^i*^ <fe. (4.10) 
I* J I 

We now decompose gj^ on the basis functions ((p J m a )o< a <d of polynomial of degree 
lower than d on each subinterval IL '■ 



(n k \ ■ = \ ~" (?i> k irfl 

with a G L 2 (I\j) solution of the system 
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Introducing the function $° m a defined on /, equal to (p 3 rn a on I° m and equal to zero 
on we have 

N 3 d 

a k ■= &' k $ J (4 II) 

m— a— 

Replacing (|4~TTj) in (|¥TTU|) leads to 

Let M"'^ be the projection matrix defined by 

(M a ^) n+hm+1 = J&n^pda. 0<,n<N t ,0<m< N r 
Then we have, for < n < Ni, 

d 



a=0 

with G#* = (G$*,..,<%*J*. 

In the special cases d = and d = 1, we obtain : 
Case d = 

In that case there is one basis function (p n = 1 on 7^ , and 

$j<p* Q d8= ( Qj i ds = (M°'°G^ fe ) n , 



with (M°'°) n+1 , ro+1 = ljiV ds, G*J„ = 4- J J3 . ds, <,n < JV*, < m < Nj. 



Case d = 1 



s-il 



In that case there are two basis functions ip l n = 1, <p t n<i \ = 2 — ^\ +1/2 on I n , and 
/,< #X,o ^ (M°'°Gf + M^Gf )„, 



4 &<i & - (M^Gf + M^Gf )„, 
with, for <, n < JVi, < m < JV,- 



(M°'°)„ +1 , m+1 = / V 1,, ds, (M^Ubm+i = 4 / 7 +V2 l f i !" +1/2 l / j ds > 

J/ " . •/ 1 k J m m k n 

(M 1,0 ) n+ i m+ i = 2 / S ~ 7 +1/2 l z , l It ds, (M 0,1 ) n+1 m+1 = 2 / 1,, S ~ J +1/2 l f< ds, 

Jl k 3 m m " Jl m K 

and G^ Q) l defined by 

"771,0 — jr ip m yj,i ab 
Aj,k _ _6_ f !Z^±1Z3 ^fc rf s 
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The projection matrices M Q '^ are computed by a simple and optimal projection al- 
gorithm without any additional grid (see [TJ, [5])- 

We now discuss the space dicretization using finite elements. 

4.2. Space discretization. We suppose that each subdomain f2i is provided 
with its own mesh 7?, 1 < i < I, such that 

For T £ T%, let h T :— su~p x yeT d(x,y) be the diameter of T and h the discretization 
parameter 

h = max hi, with hi = max hx- 

i<i<i TeT h * 

Let V\ (T) denote the space of all polynomials defined over T of total degree less than 
or equal to 1. Then, we define over each subdomain the conforming spaces V? by : 

Vl = {v i>h e C°(n,), »i, fc | T € Pi(T), VT e it). 

In what follows we assume that the mesh is designed by taking into account the 
geometry of the T^ j in the sense that, the space of traces over each j of elements 
of Vfl is a finite element space denoted by V 1 ^ . Let n lJ be the dimension of Vlf and 
(Xp''h,)i<i<n i <j the finite element basis functions of V^ J . 

We consider two cases : when the grids in space are conforming, and the case of 
nonconforming space grids. 

4.2.1. Conforming case. In the case of conforming grids in space, we have 
VV 3 = V£\ We can replace Bi by V£ in the variational formulation. We set : 

di(u,v) = I (-((bi-Vu)v— (bi-Wv)u)) dx + I ViVu-Vvdx + I (c, + — V -bijuv dx, 

(4.12) 

and 

< C it jU, v >r i , J = J r . . ((Pi,j ~ ^y 1 ) uv 

+q it j (d t (l l u) + Vrv 3 • (riju))v - s;,jVr 4(i uVr <a 'i>) da, 

< Ciju, v >r < , 3 = J" r . . ({Pi,j - h Y L ) uv 

+q it j (d t {l j u) + Vr i>3 - ■ (ri t ju))v - s 1J Vr, 3 MV r , J i') da. 

We introduce the discret algorithm : let (gj ^ h ) be a given initial guess in Vdiy]l J ,%), 
for 1 < i < 7, j € Afi- Let Uf h be the approximation of u\ in ~Pd(VZ,Tj). Then, at 
iteration k > 1, we solve the subdomain problem in Qf. 

(d t (TU£ h )vi, h + ai(Ul° >h ,Vi, h )) dx+< CijUftoVifi > r , j 

r P l fv hh dx+ [ gi jh v ijh da, in (0,T), Vv i>h e V£, (4.13) 

For k > 2, Vh € V£ J , we define 

/ 9 k i,j,hVhda := P i j g^, h v h da, (4.14) 
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r 4 



with 

9j, ilh Vh da : j g)~lv h da+ < CiiUfa 1 + C hJ U^\v h -i ■ 

In equation (|4.14p we used the fact that the space of traces over each Tij of elements 
of VI is the same as the space of traces over each Y^j of elements of V£ . For the 
computation of the right-hand side in (I4.13[) . we follow the same steps as in section 
14.1.21 where we replace with g^ i h € V]f defined by 



9j,i,h = ( / 9j,i,hXii dv, •••• / 9j,i,hX%,S h da ) > 



and we replace G l ^ a with a h £ V]f solution of 

E^!a,h /. ^k a <pds= f /9e{0,...,d}. (4.15) 



The discrete formulation in the cases d = and c? = 1 are obtained from (|4.5I) and 
(|4.6p. by replacing Bj by VjJ. 



When the space grids are nonconforming, following |8|, we cannot replace directly 
Bi by the finite element space Vi in the variational formulation. We have to consider 
equation (|4"Tj) (i.e. {0| for d = 0, and for d = 1). 

4.2.2. Nonconforming case. In this section we extend the nonconforming ap- 
proach in [SJ. We consider the mortar spaces Wl' 1 as in [8]. Let m IJ be the dimension 
of Wlf and (V'fe"/ l )i<fc<m i .j the finite element basis functions of Wl J . We introduce 
the discrete algorithm : let (U^ 1 , Q^ 1 ) € "PdO^, 71) X VdiW^ ,%) be a discrete ap- 
proximation of (U^ , VidmUi 1 ) in fij at step fc— 1. Then (C/^, is the solution 
mP d (VZ,%)xP d {Wl' j ,%) of 

^(^C*fc.«i,h)»- + 5i(C^fc»«i,/.)< 

+ ^ (Qi,h - b -^Ul h )v l , h d<7 = (P l f,v hh ) t , in (0,D, Vu i>ft e v& 

/ (Qa - ^ • n i u i,h + Pij u* h )rtf da 

+ I (^(9 t (I l ^J + V r ,, j .(r lJ ^J)^+^ 5jJ Vr i X/ l Vr l , J ^)^ (416) 



r, 



QW - b, ■ n^- 1 + Pl , U^ 1 )^ da 
on (0.7). V^ j 'g WjfJefA. 



We give first the interior scheme for d — and <i = 1 and then the computation of 
the right-hand side in the transmission condition of (|4.16[) . 
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Interior scheme. The discrete problem in subdomain O = f2j in (|4. 1 6[) is defined 
as follows : find (U h ,Q h ) := {U^Q^) in V d {V£,%) x V d (W^ ,%) solution of 

d P b Th 

- (2U h , v h ) + a(U h , v h ) + j^(Q h - —U h )v h da = {Pf, v h ), in (0, T), Vv h € 

J ({Q-b-nU h + P U h + q(d t (lU h ) + V r • {rU h )))^ h + qsWvU h V v ^h) da 
= J r (Pg)iP h da, on (0,T), G 
In the cases d = and d= 1 we obtain 
Case d = 

In that case, the approximating functions are piecewise constant in time: Uh{t) = 
U£ +1 — Uft, and Qh(t) = Qh +1 — Qh + 011 Ini an d the discrete problem reduces to 
find {U% + \Ql +1 ) e Vl x W l h ' j solution of 

rrn+1 /■ h 7i 

(^,v h ) + a(uj; + \v h ) + / (q;: +1 - — ^ +1 ) Uh d<7 
r/ n i f 

"Vi "Vi ■//„ 

/ - 6 • n[^ +1 +pf/^ +1 + + V r • (r£/£ +1 ))V/i + «aV r l# +1 Vr^) da 

Jr "Vi 

= I 9iH>hdiT+±- [ f g(s)ip h dads,W h £WZ j . 
Jr "Vi j/ n Jr 

Case d = 1 

In that case, we write U h (t) = U^ 1 + 2 *~^+ 1/2 Lffi 1 and Q h {t) = Q^ 1 + 
2 t ~^+ 1/2 Q^+ 1 on li,U% = U h (t n ), and the discrete problem reduces to find (U^\ Q^ 1 ) 
and (E/"£\ Q™^ 1 ) in VJj x V^ J solution of the system 

^(KV + U?£\v h ) + ~a(U^ h \v h ) + jf (Q~W ~ ^USX^Vh da 

"Vi "Vi J/„ 

+ U?£\w h ) + a{U#\w h ) + Jwti 1 h -^U^)w h da 

= -g,^) + A [ r s ~[ n+1/2 (f, Wh )ds, v Wh evi, 

[ ~ & • nf/^ 1 + f (C/^ 1 + U^) + qVr • (rE/ft 1 ))^ da 

Jr "'n 

+ / qsVrU^V^h dtr = / <Z^/i da + ^ [ [ gyj h da ds, V^/, € W™, 
Jr Jr "Vi fc n j/„ Jr 

+ / qsVrU^VrCh da = - [ q^( h da + A / *" +1/2 / <?0, dff *», V& € 
Jr Jr "Vi "Vi J/„ "Vi Jr 

Transmission terms. Let (U9 h ,Q° h ) e V d (V^,%) x V d (Wi' j ,71) be a given 
initial guess, for 1 < i < I. Then, at iteration k > 1, we solve the subdomain problem 



Schwarz Waveform relaxation and discontinuous Galerkin 



27 



in 

d 



hi ■ n 



(Qth ~ ^U£ h )v i<h da = (P7,w Ml k in (0,T), V« ilh G 
(Qi,h " b, ■ n,i :: h + PiJ Ul h )^ da (417) 

P i 9h((Ufa\Q k j - h 1 ),1>i i ), on(0,T), V^e^jeM, 



r, 



with, for fc > 1, 



g h ((U, Q),tp) := . (-Q + 6^ • n^C/ + C/)^ d<r 

+ f QijidtiFU) + Vr M ■ (rijU))^ + q l , 3 s l , J V T ^UV T ^) da. (418) 



For the computation of the right-hand side in (|4.17[) , we follow the same steps as in 
section l4.1.2l where we replace with 

5^ = (^((^^^ 1 )^ll) I ■•■^,((^^Q^ 1 )>;^, J ,J)^ 

and we replace G 1 ^ a with G 1 ^ h € Wt^ solution of 

E^ca /. <pL,«<pL,f>da = f fo^p, p e {0,...,d}. (4.19) 

For the computation oi g^ i h , we write (U^ 1 )^, . = X)™=i u j.iXi'.h> ana ^ Qj,^ 1 = 

^1,1 ^d introduce Q^ 1 = z^f, U^ 1 = (u^, ..,«£„,.,)', and 

the projection matrices, for 1 < fc < m'' 3 , 1 < I < to-'' 1 , and 1 < £ < n^' 1 , 

Then 

#M - -M^'Q^ 1 + d^M^ 1 ) + (M& + B« + K^)^ 1 . (4.20) 

The projection matrices are computed using the projection algorithm in [5]. 

5. Numerical Results. We have implemented the algorithm with d — 1 and 
Pi finite elements in space in each sub-domain. Time windows are used in order to 
reduce the number of iterations of the algorithm. For the free parameters defining 
Sij and Si j, we chose r^j to be the tangential component of the advection, Sij the 
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value of the diffusion in the domain Qj. The optimized parameters and Qij are 
constant along the interface. They correspond to a mean value of the parameters 
obtained by a numerical optimization of the convergence factor [5] . 

We first give an example of a multidomain solution with discontinuous variable 
diffusion, for two subdomains and one time window. The advection velocity is also 
discontinuous, taken normal to the interface in one subdomain, and tangential to the 
interface in the other subdomain. The latter case of a flow tangential to the interface 
is difficult when the interface conditions are not related to the convergence factor of 
the domain decomposition method (see for example \l'2\). 

The physical domain is f2 = (0, f ) x (0, 2), the final time is T = 1. The initial value 
is M = O^Se- 100 ^- - 55 ) 2 ^- 1 - 7 ) 2 ) and the right-hand side is / = 0. The domain 
£1 x (0, 2) is split into two subdomains fli = (0, 0.5) x (0, 2) and fl 2 = (0.5, 1 ) x (0, 2). 
The reaction c is zero, the advection and diffusion coefficients are b\ — (0, — f ), v\ = 
0.001^/y, and b 2 = (— 0.1,0), v 2 — O.lsin(xy). The mesh size over the interface and 
time step in Q,\ are hi = 1/32 and ki — 1/128, while in fl 2 , h 2 = 1/24 and k 2 — 1/94. 
In Figure lSTTI we observe, at final time T = 1, that the approximate solution computed 
using 3 iterations (right figure) is close to the variational solution computed in one 
time window on a time conforming finergrid (left figure). 



Monodomuin solution. At lime L=T=I DG-OSWR Solution. At 1lmc E=T=I 




Fig. 5.1. Variational (left) and nonconforming (right) DG-OSWR solutions 



We analyze now the precision in time. The space mesh is conforming and the 
converged solution is such that the residual is smaller than 10~ 8 . We compute a 
variational reference solution on a time grid with 4096 time steps. The nonconforming 
solutions are interpolated on the previous grid to compute the error. We start with 
a time grid with 128 time steps for the left domain and 94 time steps for the right 
domain. Thereafter the time step is divided by 2 several times. Figure 15.21 shows 
the norms of the error in L°° (I ; L 2 (f2i)) versus the number of refinements, for both 
subdomains. First we observe the order 2 in time for the nonconforming case. This 
fits the theoretical estimates, even though we have theoretical results only for Robin 
transmission conditions. Moreover, the error obtained in the nonconforming case, in 
the subdomain where the grid is finer, is nearly the same as the error obtained in the 
conforming finer case. 
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Fig. 5.2. Error between variational and DG-OSWR solutions versus the refinement in time 

The computations are done using Order 2 transmission conditions. Indeed, the 
error between the multidomain and the variational solutions decrease much faster with 
the Order 2 transmissions conditions than with the Robin transmissions conditions 
as we can see in Figure [5"UI in the conforming case. 




5 10 15 20 

Iterations 



Fig. 5.3. Convergence history for different transmission conditions 



We now consider the advection-diffusion equation with discontinuous porosity uj: 

Lodtu + V • (bu ~ vVu) = 0. 

The physical domain is O = (0,1) x (0,2), the final time is T = 1.5. is split 
into two subdomains. The interface T is parametrized with a Hermite polynomial 
(i + ((2 S -l) 3 + 2(2 S -l) 2 + (2 S -l))l s <i+((2,s-l) 3 -2(2 S -l) 2 + (2 S -l))l s >i),0< 
s < 1, see Figure [531 The advection, diffusion and porosity coefficients are 
bi = (-sin(f (y-l))cos(Tr(x-l)),3cos(l(y-l))sin(Tr(x-l))), v t = 0.003, un = 0.1, 
b 2 = h, v 2 = 0.01, u 2 = 1. 
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Fig. 5.4. Domains Qi (left) and Q2 (right) 

We first consider a conforming grid in space. The time step in fii is k± = 1/180, 
while in 0,2, &2 = 1/100. In Figure 15.51 we observe, at final time T — 1.5, the 
approximate solution computed using ten time windows and 5 iterations in each 
time window. It is close to the variational solution computed in one time win- 
dow on the conforming finer space-time grid as shown on the error, in Figure 15.61 




□ D2 



OS OS 1 




Fig. 5.5. DG-OSWR solution after 10 
time windows and 5 iterations per window 



Fig. 5.6. Error with variational solu- 
tion after 10 time windows and 5 iterations 
per window 

We analyze in Figure 15.71 the precision versus the time step. The converged 
solution is such that the residual is smaller than 10~ 12 . A variational reference solution 
is computed on a time grid with 7680 time steps. The time nonconforming solutions 
are interpolated on the previous grid to compute the error. We start with a time grid 
with 120 time steps for the left domain and 26 time steps for the right domain and 
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divide by 2 the time steps several times. Figure 15.71 shows the norms of the error in 
L°°(I; L 2 (S!j)) versus the time steps, for both subdomains. 




We observe the order 2 in time for the nonconforming case that fits the theoretical 
estimates. In Figure l5"7H we show also the norms of the error in L 2 (Q,i) at final time 
t = T versus the time steps, for both subdomains. We observe the order 3 for the time 
nonconforming case. This corresponds to the superconvergence behavior described in 

We now consider nonconforming grids in space as well. The mesh size and time 
step in are hi = 0.032 and k x = 1/120, while in ti 2 , h 2 = 0.048 and fc 2 = 1/26. In 
Figure f5T8l we observe, at final time T = 1.5, that the approximate solution computed 
using 5 iterations in one time window is close to the solution computed with the 
conformal in space grid in Figure 15.51 left. In Figure 15.91 and 15.101 we observe the 
precision versus the mesh size. The converged solution is such that the residual is 
smaller than 10~ 12 . A variational reference solution is computed on a time grid 
with 960 time steps and a space grid with mesh size h = 3.5 10~ 3 . The space-time 
nonconforming solutions are interpolated on the previous grid to compute the error. 
We start with a time grid with 60 time steps and a mesh size hi = 0.056 for the left 
domain and 20 time steps a mesh size hi =0.11 for the right domain and divide by 2 
the time step and mesh size several times. Figure [3~9l shows the norms of the error in 
L 2 (I; L 2 {D,i)) versus the time steps, for both subdomains. We observe the order 2 for 
the nonconforming space-time case, even though we have theoretical results only for 
the time semi-discrete case. Figure 15. 101 displays the norms of the error in L 2 (£li) and 
in 7? 1 (r2i) at final time t = T versus the mesh size, for both subdomains. We observe 
the order 2 for the L 2 error, and the order 1 for the H 1 error for the nonconforming 
space-time case. 
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Fig. 5.9. Relative L 2 error in time and Fig. 5.10. I? and H 1 errors at the final 

space time 

6. Conclusion. We have proposed a new numerical method to solve parabolic 
equations with discontinuous coefficients. It relies on the splitting of the time interval 
into time windows, in which a few iterations of an optimized Schwarz waveform relax- 
ation algorithm are performed by a discontinuous Galerkin method in time, with non 
conforming projection between space-time grids on the interfaces. We have shown 
theoretically in the Robin case that the method preserves the order of the discon- 
tinuous Galerkin method. Numerical estimates of the L 2 (I; L 2 (Hi)) error and the 
H 1 error at final time have shown that the method preserves the order of the space 
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nonconforming scheme as well. The analysis of the fully discrete scheme will be done 
in a further work. 
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